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INTRODUCTION 

This  report  discusses  the  design  and  implementation  of  the  adaptive 
three  dimensional  maneuvering  target  tracking  filter  in  spherical  coor- 
dinates(  developed  at  Virginia  Polytechnic  Institute  and  State  University 
over  the  past  two  years. 

v All  filter  performance  statistics  concerning  intercept  range, 
probability  of  kill,  cumulative  probability  of  kill,  S1G3D,  etc.  were 
obtained  from  the  statistical  analysis  program  called  GUSS.  This  program, 
supplied  by  the  Naval  Surface  Weapons  Center,  Dahlgren,  was  installed 
on  the  VPI&SU  computing  system  for  this  purpose. 

Since  GUSS  performs  its  analysis  in  rectangular  coordinates  while 
the  filter  filters  in  spherical  coordinates,  two  additional  sets  of 
subroutines  were  required  to  act  as  buffers  between  GUSS  and  the  filter. 

The  first  of  these  sets  converted  the  GUSS  supplied  three  dimensional 
target  position  data  from  rectangular  to  spherical  coordinates.  The 
second  set,  which  was  much  larger  than  the  first,  converted  the  filtered 
estimates  of  target  position  and  velocity  in  each  of  the  3 spherical 
coordinate  directions  into  their  equivalent  values  in  rectangular  coor- 

j 

dinates  for  subsequent  processing  by  GUSS. 

1 

Unless  otherwise  specified,  reference  in  this  report  to  a predictor 
refers  to  the  predictor  used  for  fire  control  purposes  and  not  to  the 
one  step  ahead  prediction  process  of  the  Kalman  filter  algorithm.  This 
predictor  takes  the  updated  estimate  of  the  target  states  at  time  k 
appearing  at  the  output  of  each  filter  channel  and  predicts  ahead  an 

I 

amount  TP  seconds  (determined  by  GUSS)  to  yield  an  estimate  of  the  future 


target  position.  This  predicted  position  is  then  used  by  GUSS  as  an 


"aim  point"  at  which  to  shoot  a shell.  By  comparing  the  actual  target 
position  TP  seconds  later  with  this  "aim  point",  GUSS  can  compute  miss 
distances  in  three  dimensional  space  as  well  as  the  probability  of 
destroying  the  target.  A constant  velocity  spherical  predictor  was 
used  in  this  report. 

Over  the  past  two  years,  two  versions  of  this  filter  were  used, 
the  second  being  an  improved  design  of  the  first.  The  performance  of 
both  designs  using  GUSS  output  will  be  analyzed  with  respect  to  each 
other  and  also  with  respect  to  the  GIP  rectangular  filter. 

The  values  of  all  filter  parameters  used  in  the  various  runs  will 


be  stipulated  as  also  will  be  whether  these  values  were  arbitrarily 
arrived  at  or  derived. 
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DERIVATION  OF  THE  THREE  DIMENSIONAL  LINEARIZED 
SPHERICAL  MODEL  FOR  A MANEUVERING  TARGET 

Exact  modeling  of  target  motion  in  spherical  coordinates  (r,  e,  6) 
usually  involves  the  simultaneous  solution  of  three  very  complex  coupled 
nonlinear  differential  equations.  The  solution  becomes  even  more  difficult 
when  forcing  functions  are  used  to  model  target  maneuvers. 

There  are,  nevertheless,  several  advantages  which  make  filtering 
in  spherical  coordinates  in  conjunction  with  air  target  radar  data, 
highly  desirable.  Chief  amongst  these  is  that  the  radar  data  itself 
is  already  in  spherical  coordinates  and  consequently  the  observations 
are  linear  functions  of  the  target  state  variables.  Therefore,  if  an 
approximate  linearized  spherical  model  for  the  maneuvering  target  could 
be  obtained,  then  the  need  to  use  a nonlinear  filter  such  as  the  Extended 
Kalman  filter  would  have  been  obviated. 

To  derive  the  linearized  spherical  model,  consider  a target  whose 
motion  in  rectangular  coordinates  is  described  by 


X * -ax  + Ujj  + w ' 

w ' ■ -aw  * + w (1) 

X XX 


where 

a is  a drag  coefficient 

u^  is  the  deterministic  input  in 
w^'  is  the  Singer  [L]  correlated 


the  x direction 

acceleration  process  acting  in  the  x 


direction  with  a time  constant  tc 


1 

a 
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w^  is  a white  Gaussian  random  process  acting  in  the  x direction. 

A similar  set  of  equations  exists  for  the  y and  z directions. 

The  significance  of  w^ ' in  relation  to  the  final  form  of  the  filter  to 
be  derived,  will  be  discussed  in  detail  later. 

Defining 


a3  x 

the  following  continuous  time  state  variable  model  is  obtained  for 
equation  (1) 


0 1 
0 -a 
0 0 
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Discretizing  (2)  in  time  yields 


1 A B 

0 E F 

0 0 e"at 


k+1 


u + 
*k 


*k 


(3) 


where 

A - (1  - e~aT)/a 

B - [L  + (ae-0*  - ae“aT)/(a  - a)  ] / (ota) 
C - (oT  - 1 + e“°T)/a2 
D - [T  + (aA  - aJ)/(a  - a)]/(aa) 


G - (J  - A) /(a  - a) 
J • (1  - e aT)/a 


A similar  state  variable  model  is  assumed  to  exist  for  the  y and  z 
directions. 

With  the  development  of  the  above  state  model  complete,  the  derivation 
of  the  linearized  spherical  model  is  given  next. 

If  (x,  y,  2)  represents  the  position  coordinates  of  a maneuvering 
target  in  rectangular  coordinates,  then  the  corresponding  position  of  the 
target  in  spherical  coordinates  is,  from  figure  1 


, 2 . 2 . 2.% 

(x  + y + z ) 

(4) 

sin  ^(z/r) 

(5) 

tan  *(y/x) 

(6) 

where  r,  e and  8 represent  range,  elevation  and  bearing  angles  respectively. 

The  linearized  spherical  model  is  derived  by  using  the  following 
expansion  [2] 

rk+l  = rk  + ^k(xk+l  ‘ xk)  + l^'k^k+l  ‘ yk>  + ff!k(zk+l  * Zk>  (7a) 

vi 3 % + ff(k(xk+i  • v + !?ik(yk«.  ■ v + lfik(zicfi  • v (7b) 

and 

0k+l  = Bk  + a^k(xk+l  " V + 3ylk(ykH-l  “ ykJ  + a^k^k+l  ~ V 


(7c) 
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r CHANNEL 


calculation: 


3r  1.  2 2 2.-%  , 

— • -(x  + y + z ) (2x)  - 


3x  2 


, 2 . 2 . 2. ^ r 

(x  + y + z ) 


3r  y , 3r  z 

— ■ and  — - — 

3y  r 3z  r 


t-|Ii+|lj+|l4  - jxx  + yy  _+  zi) 
3x  3y  7 3z  r 


Substituting  these  relations  into  7(a) 


rk  + r'k(xk+l  " xk}  + rik(yk+l  “ yk)  + r'k(zkfl  “ zk> 


From  (3) 


(Xk+l  - Xk>  " ^k  + BVx^  + Cuxk  + Dwxk 

(yk+l  • yk}  ■ Ayk  + BWy'  + % + (9) 

7k  7k  7k 

(z  - z ) - Az.  + Bw  ' + Cu  + Dw 

k+i  k k zk  zk  Zk 

•*  rk+l  " rk  + ?lk[Axk  + Bvx^  + Cuxk  + ^xj  + 

*|  [Ay  + Bw  1 + Cu  + Dw  ] + 

K X yk  yk  yk 

rlkfAik  + Bwzk  + Cuzk  + Dwzk1 

- rk  + A[^-±  zz]  |k  + B[vx-(^+  wy’($  + vz’(f),  |k 

+ C[ut(^+  uy(r>+  Virile  + D[wx(T)+  wy(f>+  wzd)J 

XV  z 

From  Figure  1,  — , £ and  — are  the  direction  cosines  between  the  X, 

Y and  z axes,  respectively,  and  the  r direction.  Therefore 
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wx'$)+  Wy'(J)+  Wz'&~  Wr' 


(10) 


Is  the  sum  of  the  projections  of  w ' , w ' and  w ' onto  the  r direction. 

x y 2 

This  sum  acting  in  the  r direction  can  be  replaced  by  an  equivalent 
single  term  denoted  by  wr' • In  a similar  manner  the  coefficients  of  C 
and  D are  called  uf  and  respectively.  Using  (8) 


r,  . - r,  + Ar.  + Bw  ' + C u + D w 

k+i  k k r k rR  rk 


(ID 


£. calculation: 
k+1 


. m xx  + yy  + zz 
r 


(8) 


In  equation  (8)  the  derivatives  of  r with  respect  to  x,  y and  z yield 
the  linear  terms  in  the  following  expansion  of  fj^ 


*k+i 3 \ + lflk(ik+i  “ V + ff^k+i ~ V + S^^k+i  ■ V 


af 


3r. 


From  (3) 


3r  x 3?  y 3f  2 
3x  " r ’ 3$  " r ’ 3z  r 


'Vi  - V ‘ + ^ + Auxk  + - *k 

«k+l  - V ■ EK  +EW,'  +A“y.+  G"v.  ' 


(12) 


<!fckl  - V • E2k  + ^ +*“zk+G“zk-  lk 
■'  fk+l  ’ fk  + ?lktE*k  + Pw^  + AU  ^ + ‘ *k] 

+ ?k[C>k  + \ + \ + \ - *k> 


I 
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+ 7lk[Eik  + Fwz’  + Au^+  Gw^-  zk] 

r.  + E["  -*  & ■■  z-*]  + F[v  '(-)  + v '(*)  + v ' (-)  I 

k r k x vr  y r z r k 


+ A[u  (-)  + u (*■)  + u (-)],  + G[w  (-)  + w (*■)  + v (-)], 
xr  y r zr  k xr  y r zr  k 


^xx  + yy  + zz 


From  equation  (8)  the  first  and  last  terms  on  the  right  side  of  the 
above  equality  cancel.  The  other  terms  involve  the  sum  of  projections 
of  rectangular  quantities  onto  the  radial  direction.  Treating  these  as 
single  terms  acting  in  the  r direction  we  end  up  with 


r,  . . ■ Er.  + Fw  + Au  + Gw 
k+1  k r,  r.  r, 
k k k 


(13) 


It  remains  to  obtain  a state  variable  model  for  w'.  Since  w'  is  a zero 

r r 

mean  correlated  Gaussian  process  acting  in  the  r direction,  a convenient 
continuous  time  model  is  given  by 


* 9 ft 

w » -aw'  + w 
r r r 


(14) 


where  w^_  is  a white  Gaussian  random  process  acting  in  the  r direction. 
Discretized  in  time  (14)  becomes 


i -aT  , .1,,  -aT. 

i ■ e w + —(1  - e )w 
r.  . i r,  a r, 

k+1  k k 


(15) 


Equations  (11) , (13 ) and  (15)  collectively  form  the  following  state 
model  for  the  r channel  of  the  linearized  spherical  model: 


r 
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J k+1 


i A B 
0 E F 
0 0 e 


-aT 


u + 
rk 


w 


e CHANNEL 


e,  , calculation: 
k+1 


1 9 

H 

1 . -1 

u]  - 

1 

du 

e * sin  {“) 

fa  lSin 

(1 

2 *5 
- u ) 

dx 

3e  

1 

t 

r (0)  - 

ifl 

-xz 

3x 

Z2> 

2 

r 

2 

pr 

(1  ~ h) 


. 2 , 2.^ 
where  p - (x  + y ) 


Similarly 


l£  . zZ* 
3y  2 
y pr 


3e 

3z 


2 4 

(i  - h) 


r (1)  - z(f) 


r_ 

P 


£l  - 


(16) 


and 


l£  v-  + — v + — 

3x  3y  y 3z 


3e 


3e  . _ ~xz 


x - 


vz 


• . K i 

V + — ^ z 


pr 


pr 


P 


■ — — 
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2 • 2 2. 

. o 2 - (xx  + yy)  2 (x  + y )z  - (xx  + yy)z 

• • 6 “ 2 . 2 2 A . 2 2 ' 2 

cr  (x  + y ) (x  + y + z ) 

Substituting  these  derivatives  into  7(b)  and  making  use  of  (9) 

-XZx 

I T owv 

‘k  ~k  ~k 


(I?) 


e.  ■ e.  + (-Xr)  I [Ax.  + Bw  * + Cu  + Dw  ] 

k+1  k vpr*-  ‘kL  K xu  xk  x,_ 


+ I . [Ay.  + Bw  1 + Cu  + Dw  ] 

Ar^  1 Lr  *■  ^ Ir  xr  \r  \r 


pr 


+ & ikiAik + bv  + cv - i 

k k k 

- + A[-(xx  +,  yy)z  + P.\  + „ - + ZZ£  w • 

k prz  k pr  x pr  y 

+ t vk'K 

+ ClT?  “*  - £ uy  + r u»Jk(r>k  + D|C‘f 
From  figure  2,  the  projection  of  the  unit  vector  in  the  direction 


of  increasing  e,  i ontothe  X axis  is  given  by 

, „ rZ\  /XN  -xz 

-sin  e cos  0 ” -(—)(—)  *■ 

Similarly,  the  projection  of  i onto  the  Y axis  is  given  by 


-sin  e sin  0 =*  -0^)  (£■)  = 

The  projection  of  ig  onto  the  Z axis  is  given  by 

cos  e » p/r 


(18a) 


(18b) 


(18c) 


In  view  of  equation  (18),  the  coefficient  of  B in  the  above  expansion 
(^)  w • + (- *5 ■)  w ’ + (£-)  w ’ 

' ' nr  xr  r 7 


of  e^,  namely 


pr 


pr 


is  the  sum  of  the  projections  of  wx'»  Wy'*  an<*  w2'  act^n8  ln  t^ie  direction 
of  ig  ; denoting  this  sum  as  we'  we  have  that 

(— ) w ' + (^-)  w ' + (£)  w ' 2 w ' 
pr  x pr  y r z e 


(19a) 


Figure 
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Similar  ly 


(Z«) 

pr 


u + (=&) 

x pr 


u + (£-)  u 
y r z 


u 

e 


and 


(z«)  w + (ZZ£) 


, _ „ , w + (^-)  w 

pr'  x pr  y vr  z 


- w 


Taking  note  of  equation  (17),  the  expansion  for  ek+1  can  now 


written 


Vl  : *k  + A\  + B(we’/r)k  + C(ue/r)k  + D(we/r)k 


e^^  calculation: 


In  equation  (17),  the  derivatives  with  respect  to  x,  y,  and  i 


the  only  ones  which  lead  to  linear  terms  in  the  expansion  of  e 


k+l‘ 


ek+l  : ®k  + 3xlk(xk+l  ' Xk}  + 3ylk(^k+l  " V + 3i'k(ik+l  ~ V 


3e  _ -xz  de  _ -yz  de 
3x  pvz ’ 3y  " prZt  3z  r‘ 


-xz  3e  -yz  3e 


Using  equation  (12)  we  get 


Vl  * \ + (^}lk  txk  + Fwx^  + Auxk  + Gwxk  - \] 

+ Ik  lE*k  + ^ + Auyk  + Gwyk  " K] 

+ (^lk  [Eik  + + AV  + GWZ,,  - V 

k k k 

■ \ ♦ E + + V 

+ F + <^>V + <f>»z'V% 

* * + + <f>-zV'k 

+ G [(— )w  + (•— ^“) w + (£)w  ].  /r, 
pr  x ' pr  y r zJk  k 

f2 i - (xx  -f  yy) z -I 
1 pr3 Jk 


(19b) 

(19c) 

be 

(20) 

are 
Thus , 
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From  equation  (17),  the  first  and  last  terms  on  the  right  side  of  the 
above  equality  cancel.  Hence  using  (19) 

Vl  " + F(we'/r)k  + A(ue/r)k  + G(we/r)k  < 


A discrete  time  model  similar  to  equation  (15)  is  used  to  describe  we': 

i —at  ■ . 1/,  -at.  f . , 

w ' - e w ' + — (1  - e )w  (2. 

k+1  k \ 

Equations  (20),  (21),  and  (221  collectively  form  the  following  state 
variable  model  in  the  elevation  plane: 


1 A (B/rk) 
0 E (F/rk) 

0 0 e-at 


+ G/rk„ 


6 CHANNEL 


0. ,,  calculation: 
k+1 


0 - tan-1  (^) 

X 


d . -1  i 1 du 
dx  fcan  ul  “ 1 + u-  dx 


ii L_  (z3P)  - J3L  .T  , =1 

3x  1 + y^  kx^  xz  + yz  p7 


i£ l ,U  . x 

3y  TT?  X 


2i-  o 

32 

i 30  . 30  . 

8-t— x + — y- 
3x  3y 


Substituting  these  derivatives  into  equation  7(c)  and  making  use  of  (9) 


Bk+1  " ek  + Ik  1 ^k  + Bwx^  + C\  + k] 

+ ^Ik  Uyk  + Bwy'  + Cuv  + ^y  1 
M yk  yk  yk 


- 6.  + A 
k 


+ c K^)UX  + £>uy]k/»t  + D [<^)wx  + <f)»yJk/ek 

From  figure  3,  the  projection  of  the  unit  vector  in  the  direction  of 

increasing  (3,  i^,  onto  the  X axis  is  given  by 

-sin  (3  - - ^ 

P 

Similarly,  the  projection  of  i^  onto  the  Y axis  is  given  by 

o X 

cos  8 “ — 

P 

In  view  of  equation  (25),  the  coefficient  of  B,  namely 

(-  *)w  « + (i)w  ’ = w ' 

P X V y 8 

is  the  sum  of  the  projections  of  w ' and  w ' onto  the  i_  direction. 

x y 8 

Denote  this  sum  by  w ' . 

8 


(25a) 


(25b) 


(26  a) 


Similarly 


(-  ^)u  + (— )u  = u. 

p x p y S 


(-  ^)w  + (— )w  = w 

P x V y B 


(26b) 


(26  c) 


Taking  note  of  equation  (24),  the  expansion  of  6,  . , can  now  be  written 
^ ■'  R+l 

Vi  ■ 6k + A\ + B<v/‘>)k + cvo,k  + ■xyo’k  <! 

3^^  calculation: 

In  equation  (24>,  the  derivatives  with  respect  to  x and  y are  the 
only  ones  which  lead  to  linear  terms  in  the  expansion  of  § . Thus, 

Bk+1  " Bk  + ailk(:Vl  ~ V + afV^k+l  " V 

38  ^y  36  x 

IT  p*’  3y  p7 
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Uslng  (12) 

Vi  ■ K * <?>  Ik 1 **k + + A\ + °\  ■ v 

+ <fr>L[  £yk  + F"  ' + *“  + o.  - yj 

p * K yk  yk  yk 
- Sk+  E ((~**p4  lk  + r K^*,’  + <*»«, ’V'k 

* 4 C)UX  + <f>uy'k/Pk  * G *>z  + ^“y'k^k 

- [-**  + 

1 p2 Jk 

From  equation  (24),  the  first  and  last  terms  on  the  right  side  of  the 
above  equality  cancel.  Hence  using  equation  (26) 

Vi  ' E“k  + F<”s'/p)k  * A<Vp)k  + G(Vp)k  <28) 

A discrete-time  model  similar  to  equation  (15)  is  used  to  describe 

* -at  , , 1/,  -aT. 

w ' - e w ' + —(1  - e )w 


w • 
0 ‘ 


'0 


k*l  6k  a ~k 

Equations  (27),  (28),  and  (29)  collectively  form  the  following  state 
variable  model  in  the  bearing  plane: 


S,. 


(29) 


8 J 


k+1 


1 A (B/pfc) 
0 E (F/pk) 


0 0 e 


-at 


6 

'(C/Pk)‘ 

(D/Pk) 

• 

6 

+ 

(A/pk) 

Ufi  + 
6k 

(G/pk) 

w * 

L 0 „ 

0 

J 

k 

■ J 

- 

w (30) 


where 


(x2  + y2) 


2\l/2 


Equations  (16),  (23),  and  (30)  constitute  the  range,  elevation  and 
bearing  channels  respectively,  of  the  linearized  spherical  model  for  a 
maneuvering  target.  It  is  observed  that  there  is  slight  coupling  between 
the  3 channels  insofar  as  the  elevation  channel  requires  the  range  esti- 
mate and  the  bearing  channel  the  range  and  elevation  estimate  to  compute  p . 
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It  should  be  noted  here  that  the  three  state  variable  models  (16), 
(23)  and  (30)  are  the  result  of  incorporating  the  Singer  correlated 
acceleration  process  into  the  linearized  spherical  model  disclosed  in  [2]. 
Summary  of  (16),  (23)  and  (30): 


— “1 

r 

• 

1 

r 

" 

w' 

LrJ 

k+l 

1 A B 

r 

"C‘ 

d" 

0 E F 

• 
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+ 
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„ „ -aT 

Ck 

0 0 e 
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(16) 


e 

"l 

A 

<*v 

•n 

"C/rk 

"D\ 

e 

- 

0 

E 

(F/rk) 

e 

+ 

A/rk 

u + 

e. 

G/rk 

w* 

k+l 

0 

0 

-aT 

e 

w’ 

_e^ 

k 

0 

k 

J 

(23) 


8 ~ 

1 

A 

<B/»k>- 

6 

5 

■ 

0 

E 

(F/Pk) 

• 

8 

+ 

(A/pfc) 

\ + 

(0/pk) 

5_ 

0 

0 

-aT 
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w' 

L_eJ 

0 

J 

k+l 

- 

- 

k 

- 

- 

With  the  development  of  the  spherical  model  now  complete, 
will  next  be  directed  to  the  adaptive  tracking  filter. 


(30) 


attention 
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ADAPTIVE  TRACKING  FILTER  IN  SPHERICAL  COORDINATES 

The  adaptive  portion  of  the  tracking  filter  is  concerned  with 
modeling  the  unknown  control  input  to  the  target.  However,  a brief 
case  study  is  first  required  since  the  filter  proposed  here  is  an 
outgrowth  of  that  disclosed  by  Moose/Gholson  [2]. 

In  the  tracking  filter  reported  by  Moose/Gholson  the  input  is 
viewed  as  coming  from  a set  of  N discrete  levels  u^\  i * 1,  2,  ....  N. 
The  maneuvering  target  is  then  represented  as  a "random  switch"  which 
arbitrarily  selects  the  target  input  from  among  this  set.  By  utilizing 
the  semi-Markov  [3]  properties  of  this  random  switching,  a set  of  N 
probabilities  W^,  i • 1,  2,  . . . , N is  computed,  where  is  the 
probability  that  the  current  input  being  selected  is  u^  . Next  a 
weighted  sum  of  these  inputs, 

d - l u(l)  W (31) 

i-1 

is  computed  and  this  quantity  d forms  the  "deterministic  input"  to  the 
Kalman  filter.  However,  in  order  for  this  method  to  work  effectively, 
many  levels  must  be  utilized.  If  an  input  is  chosen  which  is  not  exactly 
"matched"  to  one  of  the  discrete  levels,  a bias  develops  in  the  filter 
estimates.  This  bias  can  only  be  reduced  at  the  expense  of  increasing 
the  number  of  discrete  levels  to  such  an  extent  as  to  Insure  that  a 
given  input  to  the  target  will  not  be  distant  from  one  of  these  levels. 

In  an  attempt  to  eliminate  the  bias  problem  arising  from  mismatched 
inputs  the  correlated  Gaussian  random  process  disclosed  by  Singer  [1] 
was  incorporated  into  the  Moose/Gholson  filter,  as  disclosed  in  the 
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previous  section.  The  resulting  estimator  consisted  of  three  Kalman 
filter  algorithms  based  on  equations  (16),  (23)  and  (30)  for  the  r,  e 
and  d channels,  respectively,  with  uf,  ug  and  u^  being  replaced  by 
weighted  estimates  (31)  of  these  inputs,  namely  u^,  ue  and  u^  using 
the  semi-Markov  statistics.  Figure  4 is  a block  diagram  of  the  elevation 
channel  of  this  modified  filter.  Figure  3 shows  the  results  obtained 
from  the  r (radial)  channel  of  this  filter  on  the  following  target 
trajectory. 

Using  a sampling  period  T ■ 0.3  sec,  a drag  coefficient  a - 0.4  and 

2 

a correlated  process  with  a standard  deviation  o ” 30  ft/sec  and 

correlated  time  constant  r -10  secs,  the  modified  filter  was  exercised 

c 

using  synthetic  data  on  a target  initially  retreating  at  a radial  velocity 
of  MACH  1*5  which  subsequently  turns  around  and  approaches  at  a radial 
velocity  of  -MACH  1*5. 

The  insert  below  the  graph  in  Figure  5 clearly  shows  the  extent  of 
the  mismatch  which  exists  between  the  levels  of  the  applied  control 
sequence  (dashed  lines)  and  the  discretized  levels  used  in  the  filter 
(solid  lines). 

The  estimated  target  radial  velocity  (dashed  lines)  is  seen  to 

have  large  oscillations  about  the  true  radial  velocity  (solid  line). 

The  oscillations  result  from  the  "states"  not  being  adequately  separated 

in  the  state  space  for  consistent  convergence  of  the  probabilities  to 

the  correct  level.  This  point  will  be  expanded  on  in  the  next  section; 

suffice  it  to  say  here  that  this  inconsistency  in  the  probabilities 

causes  the  estimate  u to  be  biased  off  a substantial  amount  from  the 

r 

actual  target  input.  This  alternate  biasing  to  the  right  and  then  to 
the  left  of  the  actual  target  input  causes  the  estimates  to  alternate 
from  being  too  large  to  being  too  small. 
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In  the  following  solution  to  this  problem,  explicit  reference  is 
made  to  the  radial  channel,  for  convenience;  however,  it  should  be 
understood  that  a similar  analysis  holds  true  for  both  the  elevation 
and  bearing  channels. 

The  heart  of  the  adaptive  filter  proposed  in  this  report  is  in 
the  forming  of  the  estimate  of  the  target  states  (in  each  channel ) 
from  a weighted  sum  of  estimates  conditioned  on  the  N individual  discrete 
levels,  rather  than  by  forming  a weighted  sum  of  the  N discrete  levels 
first  and  then  computing  the  estimate.  This  difference  can  easily  be 
seen  by  comparing  Figures  4 and  7. 

To  this  end  consider  the  state  model  (16).  This  state  model  views 
the  target  input  acting  in  the  radial  direction  as  being  derived  from  a 
correlated  Gaussian  density  having  a mean  value  uf.  Next  consider  a 
series  of  N such  Gaussian  curves  with  displaced  mean  values  u 
i - 1,  2,  . . . , N and  partially  overlapping  ’’tails"  as  shown  in 
Figure  6.  If  a bank  of  N Kalman  filters  is  formed,  each  filter  based 
on  (16)  with  the  deterministic  input  u^  being  a different  one  of  these 
N mean  values,  then  a series  of  N estimates  is  obtained,  each  conditioned 
on  a different  Gaussian  curve  of  Figure  6.  Next  a weighted  sum  of  these 
estimates  is  obtained  in  a manner  to  be  disclosed  below,  and  this  weighted 
sum  is  taken  to  be  the  unconditioned  estimate  of  the  target  states. 


Calculation  of  weighting  coefficients: 

As  has  been  disclosed,  the  target  input  is  viewed  as  coming  from 
one  of  N possible  overlapping  Gaussian  curves  each  of  which  has  a predeter- 
mined mean  value.  As  the  target  executes  a series  of  evasive  maneuvers 
in  the  radial  channel,  for  example,  the  changing  input  to  produce  these 
maneuvers  is  now  viewed  as  randomly  switching  among  these  N curves.  By 
applying  the  semi-Markov  statistics  to  this  switching  process  a series  of 
N probabilities  W , i * 1,  2,  . . . , N is  generated  where 

= Pr  {target  input  is  being  derived  from  the  Gaussian 
curve  whose  mean  value  is  u 

r 

These  are  then  used  to  form  the  weighted  estimate. 

Before  deriving  the  recursive  form  for  W,  i«l,  2,  ...,N 
it  should  first  be  pointed  out  again  that  each  of  the  separate  signal  channels 
is  independent  in  the  sense  that  the  processing  of  a given  channel  depends 
only  on  past  estimates  in  a parameterized  manner.  For  instance,  the 
elevation  channel  carries  on  its  processing  of  the  elevation  measurements 
dependent  on  the  other  channels  only  to  the  extent  that  the  previous  range 
estimate  is  required  as  a parameter  in  the  elevation  channel  transition 
matrices.  Because  of  this  decoupled  property  of  the  algorithm,  the 
estimation  algorithm  for  a single  (elevation)  channel  need  only  be 
described,  the  other  channels  (range  and  bearing)  being  similar. 

We  begin  with  the  well-known  relation  that  the  optimal  estimate 
of  the  elevation  signal  can  be  written  as  a weighted  sum  of  the  input- 
conditioned  estimates.  Thus,  if  'e^(k  + 1)  represents  the  optimal 
estimate  of  e(k  + 1)  given  that  the  ith  input  force  u^  is  present 
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(the  ith-input  force  being  one  of  the  previously  described  mean  values), 
then  based  on  the  data  sequence 

Z(k  + 1)  - {z(l),  z(2) z(k),  z (k  + 1)}, 


we  define 


3(k  + 1)  - l e^(k  + 1)  W (k  + 1) 
i-1 


(32) 


where 


Wi(k  + 1)  - Pr  (u(k)  - u(i)|Z(k  + 1)} 


(33) 


and 


§(i)(k  + 1)  - E{e(k  + l)|u(k)  = u(i),  Z(k  + 1)}. 


Equation  (32)  is  a total  probability  expression  developed  from  the 
basic  relation  that 

e(k  + 1)  - E{e(k  + 1) | Z (k  + 1)} 

is  the  optimal  mean-squared  estimate.  It  is  well  known  that  the  optimal 
input-conditioned  estimates  are  provided  by  suitably  matched  Kalman 
filters.  In  particular, 

a(1)(k  + 1)  » #(k)  e(i)(k)  + T(k)  u(i)  + K(k  + 1)  [z(k  + 1) 

- H$(k)  e(i)(k)  - Hr(k)  u(i)]  (34) 


where 


M(k  + 1)  - $(k)  P(k)  4T(k)  + 'i'(k)  QYT(k) 


(35) 
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K(k  + 1)  - M(k  + 1)  HTlHM(k  + 1)  HT  + R]"1, 


and 


(36) 


(37) 


P(k  + 1)  - [I  - K(k  + 1 ) H]  M(k  + 1). 

The  matrices  ♦ , T and  V are  used  to  denote  the  respective  coefficient 
matrices  in  (23);  in  addition,  due  to  the  spherical  approximation  the 
measurement  matrix  assumes  the  simple  form  (for  each  channel) 

H - [1  0 0] . 

The  following  is  an  outline  of  the  analysis  given  in  [3]  to 
calculate  the  recursive  weighting  coefficients  W^,  i ■ 1,  2,  . . . , N. 
Defining  Z(k  + 1)  ■ (Z(k) , z(k  +1)},  apply  Bayes  Theorem  to  (33)  and 
obtain 


Wi(k  + 1)  - 


Pr  (u(k)  - u(i)|z(k)>  p{ z (k  + l)(u(k)  - u(i),  Z(k)} 

p { z (k  + 1) | Z(k) > " (38) 


The  denominator  is  independent  of  i and  is  therefore  common  to  each 
W^(k  + 1)  as  a normalizing  constant.  The  first  numerator  factor  of 
(38)  is  determined  from  the  semi-Markov  input  process.  Expanding  this 
factor  in  a total  probability  expression 

N 


Pr{u(k)  - u(i) | Z(k) } - 1 Pr(u(k)  - u.|u(k  - 1)  - u. , Z(k)}W 

j -1  1 J J 

And  since  Z(k)  has  no  influence  on  the  Markov  state  transitions, 

N 


(k) 


Pr{u(k)  - u(i)|Z(k)}  - 1 & W, 

j-1  3 J 


(k) 


(39) 


where 


9 j ^ * Pr{u(k)  - u1|u(k  - 1)  - u^) 
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Combining  (38)  and  (39) 

m N 

W (k  + 1)  - C.  p{z(k  + l)|u(k)  - u'1J,  Z(k)>  l 6 W (k)  (40) 

j-l  J 

is  the  desired  recursive  relation  for  W^.  The  required  density  p is 
approximately  normally  distributed  and  has  distribution 

P(z(k  + l)|u(k)  - u(i),  Z (k) } 'b  N {mt(k  + 1),  C^k  + 1)},  (41) 

where 

m^k  + 1)  - H[*(k)  e(i)(k)  + T(k)  u(i)(k)]  (42) 

and 

Ct(k  + 1)  - [HM(k  + 1)  HT  + R]  (43) 

The  final  estimation  algorithm  consists  of  the  calculations 
implied  by  (34),  (41),  (42),  (43),  (39),  (40)  and  (32). 


With  the  derivation  of  the  adaptive  estimator  now  complete,  several 
comments  are  £n  order. 

The  use  of  the  word  "optimal"  in  the  previous  analysis  needs  to  be 
qualified  for  the  following  two  reasons: 

(1)  Since  the  target  input  is  usually  unknown  when  using  actual 
radar  data,  the  above  modeling  of  this  unknown  input  is  at 
best  approximate  and  generally  does  not  match  the  true  target 
input  from  iteration  to  iteration.  Indeed,  being  cognizant 
of  this  mismatch,  the  Singer  process  was  incorporated  to  raise 
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the  estimator  uncertainty  concerning  the  input  and  in  this 
manner  to  produce  improved  estimates.  In  reference  [4]  it  is 
shown  that  one  additional  covariance  term  is  also  needed  to 
account  for  thia  mismatch. 

(ii)  When  a target  maneuver  occurs,  the  weighting  coefficients  do 

not  respond  immediately  but  rather  have  a finite  learning  time. 
Consequently,  during  this  learning  period  the  incorrect  filter 
is  being  weighted  the  most  while  the  filter  closest  to  the  new 
target  configuration  is  being  weighted  by  a small  amount  causing 
the  estimates  to  lag  the  true  target  states. 

Because  of  (i)  and  (ii)  the  adaptive  estimator  developed  here  must  be 
considered  sub-optimal. 

Consider  the  measurement  density  conditioned  on  the  ith  mean  value 
as  given  in  (41) . This  density  has  covariance 

Ci(k  + 1)  - [HM(k  + 1)  HT  + R]  (43) 

where  M(k  + 1)  is  given  by  (35).  What  characterizes  the  different 
target  "states"  is  the  set  of  Gaussian  curves  used  to  model  the  switching 
input.  However,  the  target  dynamics  remain  the  same  for  all  the  "states”. 
Consequently,  if  the  process  and  measurement  noise  covariances  Q(k  + 1) 
and  R(k  + 1),  respectively,  are  assumed  to  remain  constant  as  the  target 
switches  from  one  "state"  to  another,  then  none  of  the  quantities  on 
the  right  of  (43)  is  conditioned  on  i.  Under  these  conditions,  for  a 
given  value  of  (k  + 1),  C^(k  + 1)  has  the  same  value  for  all  values  of 

N.  Indeed  it  is  clear  that  for  each  value  of  (k  + 1), 


i - 1.  2,  . 


• » 
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the  entire  covariance  analysis  is  identical  for  each  filter  in  the 
previously  mentioned  filter  "bank".  Therefore,  the  bank  of  filters  may 
be  reduced  to  a single  filter  for  this  segment  of  the  algorithm. 

Consider  next  augmenting  both  the  state  predict  and  state  update 
estimate  vectors  of  this  single  filter  to  form  two  N column  matrices,  the 
ith  column  of  each  matrix  being  the  estimate  conditioned  on  the  ith 
mean  value.  If  the  scalar  measurement  residual  is  augmented  to  form 
an  N element  measurement  residual  vector,  then  this  single  "augmented" 
filter  now  produces  a sec  of  N conditioned  estimates  of  the  target 
states  as  if  an  entire  bank  of  N filters  had  been  executed.  These  N 
estimates  are  then  weighted  by  the  weighting  coefficients  (40)  to  yield 
an  unconditioned  estimate.  A block  diagram  of  this  "augmented"  filter 
is  given  in  Figure  7. 

The  underlying  causes  of  the  oscillations  in  Figure  5 can  now  be 
explained.  In  equation  (42),  at  time  (k  + 1)  the  mean  value  of  the 
measurement  conditioned  on  u^  is  seen  to  be  a function  of  the  quantities 
g^(k)  and  u^(k).  Both  of  these  quantities  have  different  values 
for  each  i and  lead  to  a good  discrimination  among  the  m^(k).  In  the 
Moose/Gholson  filter,  the  corresponding  m^(k)  has  the  value 

mJ(k)  - H[$(k)  e(k)  + T(k)  u(i)(k)]  (44) 

In  (44)  S(k)  has  the  same  value  for  each  i and  only  the  quantity  u^(k) 
is  used  to  differentiate  among  the  different  m^(k) . The  improved  filter 
performance  resulting  from  the  better  discrimination  afforded  by  (42) 
can  be  seen  from  Figure  8.  The  trajectory  in  this  Figure  is  identical 
in  every  respect  to  that  of  Figure  5.  As  can  be  seen  from  the  insert  in 
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Figure  7.  Elevation  channel  block  diagram  of  the  augmented  filter. 


Augmented  Filter 


Figure  8.  Radial  velocity  tracking  by  the  augmented  filter 
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Figure  8,  the  mean  values  used  in  Che  augmenced  filter  are  also  identical 
to  chose  used  in  Figure  5;  in  addition  process  and  measurement  noise 
covariances  were  the  same  in  both  runs.  The  oscillations  experienced  in 
Figure  5 are  seen  to  be  absent  in  Figure  8.  Maneuvers  are  detected 
sooner  and  the  filter  settling  time  is  much  shorter. 

In  the  next  section  the  design  process  used  in  assigning  values 
to  the  various  parameters  in  the  augmented  filter  will  be  discussed. 
Following  that  the  filter  performance  against  actual  air  target  radar 
data  will  be  illustrated  graphically  from  the  point  of  view  of  the  GUSS 
supplied  SIG3D  parameter. 


PARAMETER  SELECTION  PROCESS  FOR  THE  AUGMENTED  FILTER 


The  parameter  selection  process  concerns  itself  with  the  following 
parameter  set 


where 


(a,  a , N,  V 

c max’ 


a} 


a is  the  assumed  drag  coefficient 

oc  is  the  standard  deviation  of  the  correlated  process 
N is  the  number  of  levels  (mean  values) 

V C^e  assumet*  maximum  possible  speed  of  the  target  set  to 

be  tracked, 
a = 1/t 

c 

The  design  process  is  given  for  the  r channel.  The  parameter 
values  calculated  for  the  r channel  are  then  used  in  the  e and  8 channels . 
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Consider  the  state  model  (16).  A continuous  time  differential  equa- 
tion for  this  model  is  given  by 

r “ -ar  + u + w ' 
r r 

w ' - -aw  ' + w 
r r r 

Concentrating  on  the  deterministic  part  of  this  differential  equation 
for  the  moment,  and  defining 


Vr  E r 


we  have  that 


V - -oV  + u 
r r r 


In  the  steady  state,  when  - 0,  the  following  relationship  exists  between 


u and  V 
r r 


u - oV 
r r 


where  V is  the  steady  state  value  of  V . 
rss 


Therefore  a bound  on  the  maximum  target  velocity  establishes  a 

bound  on  the  maximum  value  of  u . 

r 

In  the  first  version  of  the  filter  developed  here,  a bound  of  ±1,200 
ft/sec  was  assumed  for  r (both  signs  are  used  to  account  for  either  a 
retreating  or  approaching  target).  In  addition  a value  of  a - 0.4  was 
selected  because  this  value  produced  a quick  transient  response 
indicative  of  a high  speed  air  target  dynamics. 

With  both  of  these  parameters  defined,  the  resulting  bound  on  u^_  is 


-480.0  < u < + 480.0 
— r - 


The  remaining  parameters  and  N are  then  chosen  in  such  a way  as  to 
embrace  the  continuum  of  possible  uf  values  lying  within  this  bound. 

For  example,  in  the  first  version  of  the  filter,  a set  of  N ■ 5 Gaussian 
curves  was  chosen,  with  the  following  mean  values 

-480.0,  -240.0,  0.0,  240.0,  480.0 
each  curve  having  a standard  deviation  of 

o - 30  (ft/sec2), 
c 

These  values  for  N,  a and  the  discrete  levels  were  arrived  at 

c 

empirically  by  exercising  the  filter  against  a variety  of  target  trajectories 

and  observing  the  overall  quality  of  filter  estimates. 

The  parameter  a ■ — where  t is  the  correlation  time  constant  of  the 
c 

Gaussian  process  had  the  value 

a » 0.1 

in  both  filter  versions.  This  value  was  suggested  by  Singer  [1]  for  a 
maneuvering  air  target. 

To  summarize,  the  following  set  of  parameters  was  used  in  the  first 
version 

T = 0.25  seconds  (sampling  interval) 
a ■ 0.4 

a • 30  (ft/sec2) 
c 

N - 5 (45) 

V » + 1,200  (ft/sec) 
max  — 

a * 0.1 

mean  values:  -480,  -240,  0,  240,  480 
An  identical  set  of  parameters  was  used  in  the  elevation  and  bearing 
channels.  Incidentally  the  above  values  were  also  used  in  obtaining  the 


results  of  Figure  8. 
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The  choice  of  the  parameter  subset 

<N,  ac) 

is  essentially  a tradeoff  of  one  parameter  against  the  other.  For  example, 
in  the  second  version  of  the  augmented  filter  the  number  of  means  values 

2 

was  increased  to  N ■ 17.  This  permitted  a reduction  of  a to  15  (ft/sec  ) 

c 

2 

in  the  r channel  and  to  25  (ft/ sec  ) in  the  e and  6 channels.  The  smaller 

a in  the  r channel  from  that  used  in  the  e and  8 channels  was  possible 
c 

by  recognizing  that  the  target  radial  velocity  is  generally  less  than  zero 
and  consequently  only  non-positive  mean  values  were  necessary  in  this 
channel.  The  entire  parameter  set  used  in  the  modified  augmented  filter 
is  summarized  in  (46)  for  convenience. 

T = 0.25  seconds  (sampling  interval) 


a = 0.4 

2 

a = 15  (ft/sec  ) r channel 
c 

2 

= 25  (ft/sec  ) e and  6 channels 


(46) 


N = 17 


V =■  -1,200  (ft/sec) 
max 

a = 0.1 

Means  values:  0,  -30,  -60,  -90,  -120,  . . . , -480  (r  channel) 

-480,  -420,  -360,  -300,  . . . , 480  (e  and  8 channels) 
The  above  analysis  shows  that  the  parameter  selection  process  is 
essentially  trial  and  error  and  consequently  no  claim  is  being  made  con- 
cerning optimality. 


4 
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COMPARISON  OF  THE  GIP  AND  AUGUMENTED  FILTER 
USING  ACTUAL  RADAR  TRAJECTORY  DATA 

The  remaining  part  of  this  report  is  devoted  to  a comparison  of 
the  augmented  filter  and  the  three-dimensional  rectangular  GIP  filter. 

This  comparison  takes  the  form  of  a series  of  graphs,  one  for  each 
trajectory,  each  graph  having  three  superimposed  plots.  Two  of  these 
plots  refer  to  the  augmented  filter  and  the  third  to  the  GIP  filter. 

The  plots  labeled  "5  levels"  and  "17  levels"  refer  to  the  augmented 
filters  having  the  parameter  sets  (45)  and  (46),  respectively. 

Each  plot  shows  the  variation  of  the  GUSS  supplied  SIG3D  parameter 
versus  range  for  the  associated  filter;  by  superimposing  three  plots  a 
convenient  visual  comparison  of  the  relative  filter  performance  is 
produced  for  each  trajectory  run.  Since  this  SIG3D  parameter  represents 
an  index  of  the  standard  deviation  of  the  three-dimensional  error  in  the 
aim  point  for  a shell  fired  at  the  target,  any  decrease  in  this  parameter’s 
value  is  highly  desirable. 

The  graphs  show  a consistently  smaller  value  for  SIG3D  produced  by 
the  "17  level"  filter  as  compared  to  that  produced  by  the  "5  level"  filter, 
for  all  values  of  range.  The  difference  between  both  of  these  values 
and  the  corresponding  SIG3D  produced  by  the  GIP  filter  is  quite  large, 
the  latter  at  times  being  several  orders  of  magnitude  larger  for  targets 
at  distant  ranges. 

In  an  attempt  to  evaluate  the  relative  reduction  in  the  SIG3D  para- 
meter individually  contributed  by  the  spherical  filter  and  the  spherical 
predictor,  the  following  series  of  tests  were  performed.  The  set  of 
trajectories  in  Figures  9 thru  15  were  run  again,  this  time  using  the 
rectangular  GIP  predictor  in  conjunction  with  the  spherical  filter. 
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Overall,  these  tests  indicated  that  the  relative  importance  of 
filter  and  predictor  is  trajectory  dependent.  Several  trajectories 
showed  only  a slight  deterioration  in  the  SIG3D  parameter  using  the 
rectangular  predictor  while  in  others  the  deterioration  was  quite 
significant . 

One  of  the  latter  cases — trajectory  //119 — is  shown  in  Figure  16. 
In  this  figure,  the  SIG3D  parameter  using  the  rectangular  predictor  is 
considerably  larger  than  the  corresponding  values  for  the  spherical 
predictor,  for  ranges  in  excess  of  2.0  K yards.  In  view  of  these 
results,  therefore,  it  appears  that  for  fire  control  purposes  the 
spherical  predictor  is  the  preferred  one  since  it  yields  consistently 
better  results  over  a broad  spectrum  of  trajectories. 
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SIG3D  versus  range  for  trajectory  #122. 
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CONCLUSION 

By  modeling  a maneuvering  air  target  in  spherical  coordinates,  a 
simplified  three  dimensional  "augmented"  tracking  filter  algorithm  has 
been  developed.  This  filter  when  combined  with  a spherical  constant 
velocity  fire  control  predictor  yields,  for  all  target  ranges,  both  a 
consistently  higher  probability  of  kill  together  with  a smaller  value 
for  the  SIG3D  parameter  than  those  produced  by  the  GIP  filter.  These 
results  were  obtained  on  a variety  of  air  target  trajectories. 

The  simplification  in  the  filter  was  made  possible  by  the  linear 
constant  coefficient  measurement  matrix  H which  not  only  has  obviated 
the  use  of  the  Extended  Kalman  filter,  but  also  has  permitted  to  a 
great  degree  the  decoupling  of  the  range,  elevation  and  bearing  channels. 

In  a separate  series  of  tests  on  these  trajectories  to  determine  the 
relative  merits  of  the  rectangular  GIP  predictor  versus  the  spherical 
predictor,  it  was  found  that  the  spherical  predictor  generally  yielded 
better  results  and  is  to  be  preferred. 

While  the  results  of  this  report  indicate  an  improvement  by 
increasing  from  5 to  17  the  number  of  levels  of  the  augmented  filter, 
such  an  increase  cannot  be  carried  out  indefinitely.  As  the  number  of 
levels  increase,  the  process  noise  variance  becomes  vanishingly  small 
and  the  filter  becomes  less  able  to  respond  to  a change  in  maneuver 
level.  Of  course  this  increase  also  imposes  a higher  computational 


load. 
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